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Abstract 

We present a comprehensive study of semiclassical phase-space propagation in the Wigner rep- 
resentation, emphasizing numerical applications, in particular as an initial-value representation. 
Two semiclassical approximation schemes are discussed: The propagator of the Wigner function 
based on van Vleck's approximation replaces the Liouville propagator by a quantum spot with 
an oscillatory pattern reflecting the interference between pairs of classical trajectories. Employing 
phase-space path integration instead, caustics in the quantum spot are resolved in terms of Airy 
functions. We apply both to two benchmark models of nonlinear molecular potentials, the Morse 
oscillator and the quartic double well, to test them in standard tasks such as computing auto- 
correlation functions and propagating coherent states. The performance of semiclassical Wigner 
propagation is very good even in the presence of marked quantum effects, e.g., in coherent tunneling 
and in propagating Schrodinger cat states, and of classical chaos in four-dimensional phase space. 
We suggest options for an effective numerical implementation of our method and for integrating it 
in Monte-Carlo-Metropolis algorithms suitable for high-dimensional systems. 

PACS numbers: 03.65.Sq, 31.15.Gy, 31.15.Kb 
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I. INTRODUCTION 



Molecular dynamics, by the spatial and temporal scales it involves, straddles the border- 
line between quantum and classical behavior. Quantum effects — above all the very concept 
of chemical reaction — are obviously crucial. Notwithstanding, a good deal of a molecule's 
external and even internal motion can be understood on purely classical grounds: precisely 
the situation for which semiclassical methods have been conceived and are optimally suited. 
Moreover, recent developments in chemical physics, experimental as well as theoretical, are 
advancing vigorously in the direction of complex time evolution and high excitations. This 
faces adiabatic techniques like even time- dependent density-functional theory (TDDFT) [1] 
and time-dependent Hartree [2, 3] with formidable challenges while it favors semiclassical 
approaches which do not suffer from limitations in this respect. 

Therefore, the modest renaissance semiclassics in the time domain (for a comprehensive 
review see [4]) is presently enjoying comes as no surprise. It includes approaches in configu- 
ration space as well as less established phase-space methods. The well-known shortcomings 
plagueging semiclassics in coordinate representation (based on WKB approximations [5] and 
the van-Vleck-Gutzwiller propagator [6-8]), such as divergences at caustics and the root- 
search problem, have largely been overcome by sophisticated refinements of the original 
method [9]. Remarkably, most of them already switch internally to mixed representations 
(combining position with momentum) or full phase space in order to smooth out divergences 
and to reduce the set of contributing classical solutions to the neighbourhood of a single tra- 
jectory, permitting the construction of initial- value representations (IVRs): uniform approxi- 
mations [10-12], semiclassical IVRs [9, 13-15] (for reviews see [16-18]), Gaussian-wavepacket 
propagation [19] with its numerous ramifications, notably the Herman-Kluk propagator [20- 
22] and Heller's cellular dynamics [23]. A similar but independent approach, the dephasing 
representation [24-26] employs the shadowing lemma for chaotic dynamics to represent the 
classical flow through a Planck cell by a single trajectory in a semiclassical approximation 
for correlation functions. 

These methods have gained widespread acceptance in practical applications but tend to be 
technically cumbersome. Their implementation in ah initio molecular-dynamics simulation 
with "on the fly" update of the electronic evolution [27, 28] has been proposed [29]. It has 
to cope, though, with the high computational costs of calculating Hessians for the potential- 
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energy surface of the electron sector — so much so that recourse is sought even at a complete 
suppression of determinantal prefactors [30]. 

Semiclassical approximations working directly in phase space promise a fresh, structurally 
more transparent approach. They are inherently free of the above problems and thus offer a 
number of tempting advantages: (i) phase-space propagation provides a natural exact IVR 
by construction, without Gaussian smoothing, (ii) it can be formulated exclusively in terms 
of canonically invariant quantities (no projection required) that (iii) allow for simple geomet- 
rical interpretations, (iv) as far as determinantal prefactors arise, they are also canonically 
invariant, and (v) being based on the density operator, not on wavefunctions, the extension 
to non-unitary time evolution is immediate, opening access to decoherence and dissipation 
[31]. 

Yet there is a price to be paid. The two representations mainly considered as candi- 
dates for semiclassical phase-space propagation both come with their specific virtues and 
vices: The Husimi function [32] is based on coherent states [33], hence closely related to 
semiclassical IVRs and Gaussian-wavepacket propagation. Coherent states possess a natu- 
ral interpretation in terms of classical probabilities but are overcomplete, contain arbitrary 
parameters, and, in order to work optimally, require to complexify phase space [22, 34]. The 
Wigner function [35-37], by contrast, provides a parameter-free one-to-one representation 
of the full density operator but encodes information on quantum coherences as small-scale 
( "sub-Planckian" ) oscillatory fringes [38, 39]. They prevent a probability interpretation and 
pose serious problems for their propagation, pointed out more than 30 years ago by Heller 
[40]. 

Owing to these seemingly prohibitive difficulties, the semiclassical propagation of Wigner 
functions still awaits being explored in its full potentiality. For a long time, it has been 
almost synonymous to a mere classical propagation of phase-space distributions, discarding 
quantum effects in the time evolution altogether [41, 42]. As a first attempt to improve on 
this so-called classical Wigner model, it suggests itself to include higher-order terms in h 
in the quantum analogue of the Poisson bracket, the Moyal bracket [38], upon integrating 
the evolution equation for the Wigner function. Such additive quantum corrections give 
rise to accordingly modified "quantum trajectories" [43-47]. However, they tend to become 
unstable even for short propagation time [48] and suffer from other practical and fundamental 
problems [49]: As we shall point out below, the very concept of propagating along single 
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deterministic trajectories is inadequate even in the semiclassical regime. 

A more convincing alternative has emerged in the context of quantum chaos where semi- 
classical Wigner functions [38, 50, 51] are appreciated as valuable tools, e.g., in the study of 
scars in wavefunctions and of spectral features related to them [52-54]. By applying semi- 
classical approximations directly to the finite-time propagator, expanding the phase instead 
of the underlying evolution equation, an appropriate phase-space propagation method for 
semiclassical Wigner functions has been achieved [55]. It features, in a geometrically appeal- 
ing manner, pairs of (real) classical trajectories and the symplectic area enclosed between 
them as their classical underpinning — not surprisingly in view of the well-known pairing of 
paths and diagrams, one forward, one backward in time, as it occurs in the propagation of 
quantities bilinear in the state vector, particularly the density operator. 

While this technique is tailor-made to time-evolve the specific class of semiclassical Wigner 
functions [38, 50, 51, 54], semiclassical approximations to the propagator proper, indepen- 
dent of the nature of initial and final states, have been constructed [56] on basis of phase- 
space path integrals [57] and of the Weyl-transformed van Vleck propagator [52]. This 
approach opens the door towards a far broader range of applications. It has already proven 
fruitful in an analysis of the spectral form factor of classically chaotic systems in terms of 
phase-space manifolds [58]. 

Here, by difference, we pretend to demonstrate its viability as a practical instrument for 
the propagation of molecular systems, including in particular numerical simulations. We 
therefore focus on its performance in benchmark models established in chemical physics 
like the Morse oscillator and the quartic double well and study standard tasks such as 
propagating Gaussians and computing autocorrelation functions. We are confident that the 
remarkable accuracy and robustness of the scheme and its unexpected numerical stability 
and efficiency suggest it as a new competitive option in semiclassical propagation. 

In Sec. II, we introduce the basic building blocks of our method, specify our two main 
approaches, via the van- Vleck propagator and via phase-space path integration, and discuss 
general properties of the semiclassical Wigner propagator, such as its geometrical structure 
and its asymptotics towards the classical limit and towards weak anharmonicity. This section 
partially corresponds to a more detailed account of material already published in Ref. [56]. 
Readers mainly interested in numerical results might proceed directly to Sec. Ill, dedicated 
to a broad survey of the performance of semiclassical Wigner propagation for prototypical 
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models of nonlinear molecular potentials, in the computation of dynamical functions like 
autocorrelations as well as of detailed structures in phase space. We also test our method 
in particularly challenging situations, involving strong quantum effects, such as the repro- 
duction of coherent tunneling and the propagation of Schrodinger-cat states, or a classically 
chaotic dynamics. Strategies to optimize the implementation of this method in numerical, 
including Monte-Carlo-like, algorithms are suggested in Sec. IV. We conclude in Sec. V with 
an outlook to immediate extensions of our work. 



II. CONSTRUCTING THE SEMICLASSICAL WIGNER PROPAGATOR 



A. Definitions and basic relations 



Alternative to its definition as a Weyl-ordered transform of the density operator [36], 
the Wigner function can be introduced as the expectation value of phase-space reflection 
operators [51] or as expansion coefficient function in a basis of phase-space displacement 
operators [59] , besides other options. For our purposes, the standard definition suffices, 
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where p denotes the density operator and r = (p, q) is a vector in 2/-dimensional phase 
space (we adopt this ordering throughout the paper but assign q to the horizontal and p to 
the vertical axis in phase-space plots). 

We shall restrict ourselves to unitary time evolution within the dynamical group U (t) = 
exp(—iHt/h) generated by a time- independent Hamiltonian H. The extension to time- 
dependent Hamiltonians is immediate but will be suppressed to avoid technicalities. The 
von-Neumann equation for the density operator, dp/dt = (— i/h)[H, p], translates into an 
equation of motion for the Wigner function, (d/dt)W(r, t) = {H w (r), W(r, £)}Mo y ai, involv- 
ing the Weyl symbol if\y( r ) of the Hamiltonian H. For h — > 0, the Moyal bracket {• , ■ }Moyai 
[36, 37] converges to the Poisson bracket. 

The time evolution of the Wigner function over a finite time can be expressed as an 
integral kernel, the Wigner propagator G w (r", r', t), 

W{r",t) = J d 2f r'G w {r",r',t)W{r',0). (2) 
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It forms a one-parameter group, implying in particular the initial condition Gy/(r", r', 0) = 
5(r" — r') and the composition law 

G w (r", r', t) = J d 2f r G w (r", r, t - t')G w (r, r', t'). (3) 

An important quantity that must not be confused with the Wigner propagator is the 
Weyl transform of the evolution operator U(t), referred to as the Weyl propagator for short 

[51], 
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It cannot be used to propagate Wigner functions as it stands, but is related to the Wigner 
propagator by a self-convolution [55-57], 



G w (r",r',t) = — -L_ J d 2 fRe^"-^ An U^,t)U w (v + ,t), 



(5) 



with f± = (r' + r" ± R)/2. It serves as a starting point for the construction of semiclassical 
approximations, based on the van-Vleck propagator (Sec. II B) as well as on phase-space 
path integration (Sec. II C). 



B. van Vleck approximation for the Wigner propagator 

A straightforward route towards a semiclassical Wigner propagator is replacing the Weyl 
propagator in Eq. (5) by the Weyl transform of the van Vleck propagator [51, 52, 57]. 
Transformed from the energy to the time domain, it reads [51, 57], 

Mr|1) = ! /r^tM. (6) 

Y V|det(Mj(r,«) + l)| 

The sum runs over all classical trajectories j connecting phase-space points r'j to r'' in time t 
such that r = f j = (r^ + rJ)/2 (midpoint rule). Mj and fij are its stability matrix and Maslov 
index, respectively. The action Sj(rj,t) = Aj(vj,t) — Hj(r,t)t, with Hj(r,t) = H^(rj,t), 
the Weyl Hamiltonian evaluated on the trajectory j (to be distinguished from H^(r,t)) 
and Aj, the symplectic area enclosed between the trajectory and the straight line (chord) 
connecting r'- to r'- [52] (chord rule, vertically hashed areas Aj ± in Fig. 1). 
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FIG. 1: Symplectic areas entering the van-Vleck-based semiclassical Wigner propagator, Eqs. (10, 
11). The vertically hashed areas correspond to the phases Aj ± of the Weyl propagators (6) accord- 
ing to the chord rule. The symplectic area (slanted hatching) enclosed between the two classical 
trajectories Tj±(t) and the two transverse vectors r'j + — r£_ and t'- + — r'-_ determines the phase 
(11) of the propagator. The classical trajectory r cl (r',t) (dashed) is to be distinguished from the 
propagation path fj(r',i) (bold) connecting the initial argument r' of the propagator to the final 
one, r". 

Substituting Eq. (6) in (5), one arrives at 

Gw( r ", r ',i) = ^E/ d2/ «^ ( '"" r>R 

3-4+ 

exp [i (S j+ (r j+ ,t) - gj_(fj_,t)) - i(fi j+ - ^Jf] 

|det[M i _(r,t) + l]det[M i+ (r,t) + l]| 1 /2 ' {) 

where indices j± refer to classical trajectories contributing to the Weyl propagators U-w(r±, t) 
in Eq. (5). The principal challenge is now evaluating the i?-integration. As it stands, 
Eq. (7) couples the two classical trajectories Tj_, r J+ , to one another only indirectly via 
R, the separation of their respective midpoints. This changes as soon as an integration 
by stationary-phase approximation is attempted, consistent with the use of the van Vleck 
propagator. Stationary points are specified implicitly by the condition r" — r' = (r"_ — + 
r'j + — r'- + )/2. Combined with the midpoint rule r' + r" ± R = r'- ± + r" ± , this implies 

r' = fJ = (iJ._ + r} + )/2 > r" = fj = (rj_ + r? + )/2. (8) 

Equation (8) constitutes a simple geometrical rule for semiclassical Wigner propagation 
[51]: It is based on pairs of classical trajectories j + , j_, that need not coincide with one 
another nor with the trajectories passing through r' and r" but must have r' midway between 
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their respective initial points r'-^ and likewise for r" . This condition is trivially fulfilled for 
identical pairs ij_(t) = r 3 - + (t), as a type of diagonal approximation. By including also non- 
diagonal terms, i.e., non-identical trajectory pairs, we add classical information contained 
in third-order terms in the action and partially compensate for the general limitations of 
the van Vleck approach on the level of wavepackets. This becomes particularly evident in 
the case of pure initial states p' = Applying the van Vleck propagator in position 

representation separately to ket and bra, allowing for off-diagonal terms in the ensuing 
double sum over classical orbits, and then transforming back to the final Wigner function 
leads to the same result as the present derivation, albeit through a more tedious and less 
transparent route. 

To complete the Fresnel integral over R, we note that 

* \S ff t) S (v t)]- j ( Mj -~* M j+ -I\_ J(M i _-M j+ ) 



<9R 2 j-v -wj 2 V + 1 M j+ + \J (Mj_ + l)(M i+ + i; 

where J denotes the 2f x 2f symplectic unit matrix [52]. Combined with the determinantal 
prefactors inherited from the van Vleck propagator, this produces 

vv // / 4^^2cos(iS7 v (r",r',t)-/f) 
G^(r ,r,t) = jjX J |det(M j+ -M 3 -_)|V2 ' (10) 

our main result for the semiclassical Wigner propagator in van Vleck approximation. The 
phase is determined by 



SJ v (v", r', t) = (v j+ - Tj_) A (r" - r') + % - Sj_ 
= f ds [1,(8) A R,-( a ) - H j+ (v j+ ) + (11) 

with fj(s) = (rj_(s) + Tj + (s))/2 and Rj(s) = rj + (s) — Tj_(s). Besides the two Hamiltonian 
terms it includes the symplectic area enclosed between the two trajectory sections and the 
vectors — r'j_ and r" + — r"_ (Fig. 1). 

In the following we list a number of general features of Eqs. (10,11): 

i. Equation (10) replaces the Liouville propagator, 

G^(r",r',t) = 6[r"-r c \r',t)}, (12) 

localized on the classical trajectory r cl (r',t) initiated in r', by a "quantum spot", a 
smooth distribution peaked at the support of the classical propagator but spreading 

8 




FIG. 2: Classical skeleton of the semiclassical Wigner propagator according to Eqs. (10,11). A set 
of initial points t'- ± (panel a, left target pattern) surrounding the initial point r' of the propagation 
path is propagated along classical trajectories (dotted lines), giving rise to a corresponding set of 
final points r" ± (right deformed target pattern). The manifold formed by their midpoints (cone-like 
structure, panel b) constitutes the "illuminated area" , the support of the semiclassical propagator. 
Inside this area, two trajectory pairs contribute to each point, one of them elliptic ("upper" shell 
of the cone, blue), the other hyperbolic ("lower" shell, red). Its boundary gives rise to caustics in 
phase space where the two trajectory pairs collapse into one. Panel c depicts the generic structure of 
the underlying phase as a function of the integration variable R in Eq. (7), an odd cubic polynomial 
of the phase-space coordinates with two extrema, a minimum and a maximum (blue symbols) 
and two saddles (red x symbols), corresponding to the elliptic and the hyperbolic trajectory pairs, 
respectively. 

into the adjacent phase space and structured by an oscillatory pattern that results 
from the interference of the trajectories involved. 

ii. The propagator (10) does involve determinantal prefactors. However, they do not 
result from any projection onto a subspace like q or p and are invariant [52] under 
linear canonical (affine) transformations [60]. 

iii. It deviates from the Liouville propagator only if the potential is anharmonic. For 
a purely harmonic potential, the two operations, propagation in time and forming 
midpoints between trajectories, commute, so that all midpoint paths fj(t) = [rj-(t) + 
ij + (t)]/2 coincide with each other and with the classical trajectory r cl (r',t). This 
singularity restores the classical delta function on r cl (r',t), see Sec. II D 3. 
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iv. By contrast to position-space semiclassics [8], the number of trajectory pairs contribut- 
ing to the summation in Eq. (10) ranges between 0, 1, and 2 following a universal 
pattern (Fig. 2): Within an "illuminated sector with its tip on the classical 
trajectory, the sum contains two trajectory pairs (four stationary points, Fig. 2c). 
In the "shadow region" outside this sector, stationarity cannot be fulfilled by real 
trajectories. Along the border there is exactly one solution (two stationary points). 
Unexpected in phase space, this caustic arises as the projection onto phase space of 
the manifold of midpoints r" = (r"_ + r'' + )/2 (Fig. 2b). 

v. The set of trajectory pairs to enter the calculation of the propagator is only restricted 
by the midpoint rule (8). However, this does not constitute a double-sided boundary 
condition since every pair fulfilling Eq. (8) for the initial points contributes a valid 
data point to the propagator: no root-search. The freedom in the choice of trajectory 
pairs can be exploited to optimize algorithms, see Sec. IV A 3. 

vi. Of the two trajectory pairs contributing to the illuminated area, exactly one is elliptic, 
associated to a pair of opposite extrema of the action, one is hyperbolic and associated 
to a pair of saddles (Fig. 2c). They can be distinguished by monitoring sgndet (M J+ — 
Mj_), positive for the former and negative for the latter. The two sets form separate 
sheets of the action (11). Within each of them, amplitude and phase of the propagator 
are slowly varying functions of r", facilitating their numerical treatment, see Sec. IV. 

vii. The propagator's oscillatory pattern encodes information on quantum coherences and 
allows us to propagate the "sub-Planckian oscillations" that characterize the Wigner 
function, solving the problem of the "dangerous cross terms" pointed out by Heller 
[40]. See Sec. HIE. 

viii. The principle of propagation by trajectory pairs is consistent with the properties of a 
dynamical group. It translates the concatenation of propagators into a continuation 
of trajectories, if the convolution (3) is evaluated by stationary phase as well. 

ix. Equations (10, 11) fail if the stationary points approach each other too closely. This is 
the case near the central peak of the propagator on the classical trajectory and along 
the caustics mentioned above. The same problem arises in the limits t — > 0, % — > 0, 
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and of weak anharmonicity. It can be overcome by means of a uniform approximation 
to the i?-integration in Eq. (7), see Sec. II D 2. 

x. As a consequence of the underlying van Vleck approximation, the propagator (10) is 
not properly normalized. Specifically, the slow decay of its oscillatory tail renders it 
nonintegrable with respect to its initial or final arguments, see Sec. Ill C 1 for quantita- 
tive details. The problem is solved by the same uniform approximation as announced 
above. 



C. Path-integral approach 

The quality of the van Vleck approximation (10) to the Wigner propagator is limited 
by its principal ingredient, the van Vleck propagator itself. Just as in configuration-space 
propagation, a comprehensive solution is provided on the basis of path integrals [61, 62]. 
In phase space an analogous theoretical framework is available since the seminal work of 
Marinov [57], largely following Feynman: Time is discretized into N equidistant sections 
t n — t' + nAt, where At = T/N and T = t" — t' is the time span to be propagated over. 
At each t n , two short-time propagators, composed of the corresponding Weyl propagators 
according to Eq. (5), are concatenated by means of Eq. (3). In the continuous limit N — > oo, 
At —7- 0, the full propagator is obtained as a double path integral comprising two distinct 
path variables r and R, inherited, respectively, from Eqs. (3) and (5), 

G w (r",f;r',0 = h~ f J Br J D#e^ s « r >>W) (13) 

with an action integral 



S{{t},{R}) = J* ds|r( S ) AR(s) + #w r(s) + ^R(s) - H w r(s) - ^R(s) j. (14) 



Only r(t) is subject to boundary conditions r(t') = r' and r(t") = r" while R(t) is free and 
may be associated to quantum fluctuations. 

For this subsection, we restrict ourselves to a single degree of freedom, r = (p, q), and to 
standard Hamiltonians H(p, q, t) = T(p) + V(q, t), T(p) = p 2 /2m, but admit an explicit time 
dependence and an arbitrary degree of nonlinearity of the potential. In this case, the Weyl 
Hamiltonian, Hw(p, q), is obtained by replacing the operators p, q in H with corresponding 
real variables, H w (p, q) = T{p) + V(q, t). 
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Returning to discrete time, we are faced with evaluating the action 
Sj V ({r},{R}) = 5^jAr n AR n + #w (r n + ^Cj ~ #w (r n - ^iCj At J 



(15) 



denoting Ar n = r n — r„_i and r n = (r„_i + r n )/2. To obtain a semiclassical approximation, 
we expand the action (15) in r n around a classical trajectory r^ 1 defined implicitly in discrete 
time by Ap d = -V'(q£,t n )At, Aq cl = (p d /m)At, and around R n = 0, n = 1, . . . , N. By 
expanding around a single trajectory, we are sampling less classical information than in 
the van Vleck approximation based on trajectory pairs. This renders our path-integral 
approach inferior in some respects detailed below. Note that the action (15) is not only 
an odd polynomial in R n , it is even lacking a linear term, hence 0(R 3 ), since according to 
Hamilton's equations of motion, the first term ~ (dr cl /dt) A R(t) = R(t) A J dH w (r d )/dr d 
cancels exactly the leading order of a Taylor expansion of H w (r cl + R/2) — H w (r d — R/2) 
with respect to R. 

Expanding to third order in all variables (there is no 0(i? 4 )-term in ^({r}, {R}) any- 
way), this leaves us with 

N 

S N ({r}, {R}) = S N ({r% {0}) + £ [ (Ag n - T"(p£)p n At) P n 

n=l 

+ (A Pn + V\t, t n )q n At) Q n + ^V'"{t, tn)Ql] ■ (16) 

The N integrations over the R n can now be performed, much like the p-integrations in the 
derivation of the Feynman path integral [61, 62], resulting in one-step propagators 

G w (r„, t n , r n _!, t n _0 = 5 [Aq n - T> n „ x At] u'^Ai [v~ 1/3 (Ap n + K'^-i At)] , (17) 

where we use shorthands T n = T(p d ) and V n = V(q d ,t n ) and define v n = Ti 2 V^' 'At/8. 

— 1/3 

The factor v n is analogous to the determinantal prefactors occurring in position-space 
path integration in that it involves higher derivatives of the potential. Similarly as in the 
van-Vleck-based Eq. (10), however, it constitutes a mathematically more benign canonical 
invariant. 

Taking into account that the one-step propagators (17) are already written in a frame 
moving with the classical trajectory, we observe three independent actions of the propagator 
on the evolving quantum spot: (i), a rigid shift along the classical trajectory, (ii), a rotation 
(shear) with the linearized phase-space flow around the classical trajectory, if it is elliptic 
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(hyperbolic), and (iii), quantum Airy- function spreading in the negative p-direction if V"'(q) 
is positive and vice versa. Their concatenation is most conveniently performed in Fourier 
space, where, under certain conditions fulfilled here [63], the convolution (3) reduces to a 



mere multiplication. With rescaled coordinates of common dimension a/ action, 

P = = {(V"/T>r 1/4 P, (V"/T")^q), (18) 
the transformation of the Wigner propagator to Fourier phase space 7 = (a, 0) is given as 



G w (Y,t'W,t') = dp" 2 J dp' 2 e-^"^'^G w (p",t";p',t'). 



Going to the continuum limit, we thus obtain 

GW', V, = - M(t", t> ) 7 '] exp 



dta(t) {M(t,t')j'l 



(19) 



(20) 



Here, M(t,t') = dr (t)/dr (f) is the stability matrix of r (t), the subscript /3 selects the 



second component of 7, and 



a(t) 



h 2 



V"[q c \t),t] 



-3/4 



V"'(q c \t),t). 



(21) 



T"[P d (t)] 

Inverting the Fourier transform (19) we recover the Wigner propagator. In the following 
subsection, we work out this result in detail for the specific case a(t) = const. 



D. Weak anharmonicity, short-time and classical limits 



1. Path-integral approach for weak cubic nonlinearity 



Equation (20) permits an analytical evaluation in the case of constant second- and third- 
order derivatives of the potential, which we will discuss in this subsection. It can arise in 
different scenarios: To begin with, consider a static binding potential with cubic nonlinearity, 
e.g., V(q) = {muj 2 /2)q 2 {l + g/3g min ), and a trajectory close to the quadratic minimum at 
9minj so that the linearized flow is elliptic. Then d 3 V(q)/dq 3 = V" = moj 2 /g min = const 

const. In accordance with Eq. (21), denote 



" /T"W2 



CO 



and we can also assume (V"/T 
a = (h 2 /8)(mui)~ 3 / 2 V'" = {ft 2 / '8q m i n ) \J '00 j 'm. In this case, the stability matrix corresponds 
to a rotation by r — u(t — t'), 



M(t,t') 



cos r — sin r 
sin r cos r 



(22) 
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and the phase of the Fourier transformed propagator (20) reduces to elementary integrals, 

_• nr" 

/ dr (a' sin r + f3' cos r) 3 . 

In a reference frame p = (fj, £) rotating by r/2 with respect to p, fj — r\ cos(r/2) — £ sin(r/2), 
£ = ?7sin(r/2) + £cos(r/2), the propagator after inverting the Fourier transform takes a 
particularly transparent form 

G^(pV; p',f ) = 7rV75)«(«)~ 2/3 Ai[«(s)- 1/3 p-]Ai[K( S )- 1 / 3 p + ], (23) 

where we abbreviate k(s) = cxs 3 /(s)/3u;, /(s) = 3s~ 2 — 1, s — sin(r/2), and p± = 
(fj ± y/ f(s) £J /2. The corresponding result for the hyperbolic case V" /T" < 0, e.g., 
close to a quadratic maximum of a cubic potential, is obtained from Eq. (23) by replac- 
ing f(s) = 3s~ 2 + 1 and s = sinh(r/2). 

Equation (23) provides us with a precise geometric outline that complements the features 
of the Wigner propagator pointed out in the sequel of Eq. (11): 

i. The quantum spot formed by the propagator fills a sector in phase space of opening 
angle 6 = 4arccot \J f(s) (Fig. 3c), rotating by r/2 if r is the angle coordinate along 
the classical trajectory. Its nodes form straight lines running parallel to its sidelines, 
with separations given by Airy-function zeros (Fig. 3). 

ii. For short scaled time r <C 1, the oscillatory tail is quasi-one-dimensional, ^ < 1, 
pointing in the negative p-direction if V" > and vice versa, and scales with time as 
t 3 . It is periodic in r with period 2ir and invariant under time reversal, r — > tq — r, 
p —> —p, with respect to r = 0, ir. 

iii. As expected for a uniform approximation, Eq. (23) resolves the sharp caustics along 
the outlines of the quantum spot present in the van Vleck approximation (10) (Fig. 3b) 
into a smooth penumbra (Fig. 3c, d). This restores the correct normalization of the 
propagator. On the other hand, it restricts the oscillatory pattern to straight node- 
lines while with Eq. (10), nonlinear deformations owing to higher nonlinearities in the 
potential can well be reproduced, see Sec. Ill B. 

iv. The finite weight the propagator (23) achieves in the penumbra zone outside the illumi- 
nated region of Eq. (10) can be interpreted as the contribution of complex trajectories. 
However, we shall not pursue this issue here. 



g w ( 7 ", V, = 6 b" ~ M (t"> t'h'] ex P 
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2. Perturbation theory from van Vleck approach 



A uniform approximation equivalent to Eq. (23) can be obtained through an alternative 
route which is particularly helpful in order to analyze the asymptotics at weak anharmonicity 
and at short time and the classical limit for the Wigner propagator: We treat the weak cubic 
nonlinearity explicitly as a perturbation of an otherwise quadratic Hamiltonian 



H(p, q) 



muo 



-q + eq 



(24) 



2m 2 

to obtain expressions for the various ingredients of the van Vleck propagator (10) through 
a perturbation expansion in e. 

From the perturbed orbits r"(r',t), we obtain pairs of trajectories r"±(t) = r' e / (r / ±f'/2,t) 
with initial points displaced by ±r'/2 from r', the corresponding centers r"(r', f', t) = (r", + 
r"_)/2, and the deviations Ar"(r',f',t) = r"(r',f',t) — f"(r',0,t) of these final midpoints 
from the endpoint r"(r',t) of the perturbed trajectory starting in r', 

(25) 



A P :'(r' ! r',t) = -^f' t F p r' + 0(e 2 ,t 4 ), 

oOJ 

A q ':(v',r',t) = J^r' t F r ,r' + 0(e\t 5 ] 



with time-dependent coefficient matrices 



2j± 
3 



rl 

3 



(26) 



Arising from a perturbation expansion of the equations of motion, Eqs. (25, 26) are reliable 
only up to leading order in r, hence their validity is also restricted in time to r < 2n. 
They constitute the classical skeleton of the quantum spot and provide the input for its 
determinantal prefactor and action. 

It is, however, more convenient to express these quantities in terms of Ar' e ', the final 
argument of the propagator relative to the classical trajectory starting in r', instead of the 
initial displacement F. By means of a few elementary phase-space transformations that 
diagonalize simultaneously the two coefficient matrices (26), the pair of quadratic equations 
(25) can be resolved for F as a function of Ar' e '. We thus obtain the propagator in the form 

2ixhf3 



[2±Ap» 2 - 12Ap"Aq" + 12f Ag //2 ] 1/4 
{cos [(3(Ar'i - Ar'Lf' 2 ] + sin [P{Ar'{ + Ar'Lf 2 } } , 



(27) 
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with $ = 2m z ' A /^ 2 %e l IH^ A and Ar^ = (1 =f V?>)^ft/mkp" ± 2y/3m/tAq". We identify 
the cosine term in Eq. (27) as the contribution of the pair of hyperbolic trajectories (saddles 
of the action, Fig. 3e) and the sine term as the elliptic contribution (extrema, Fig. 3f), cf. 
remark (vi) after Eq. (11) and Fig. 2c. 

Unfortunately, Eq. (27) becomes spurious precisely in the limit e — > for which it has been 
devised, since in this limit, the stationary points of the action coalesce and the stationary- 
phase approximation underlying Eq. (27) breaks down. The problem can be fixed, however, 
in a straightforward manner: Since the action in Eq. (7) is odd in the integration variable 
R and the two known trajectory pairs entering Eq. (27) mark its four stationary points 

R = — R++ and R-+ = — R+ (Fig. 2c), corresponding resp. to the elliptic and hyperbolic 

pairs, we can uniquely reconstruct an effective action that contains only linear and cubic 
terms in R (Fig. 2c) and thus permits an exact evaluation of the /^-integration. The result 
(cf. Fig. 3c), replacing Eq. (27), reads 

^Ai (^) Ai (*± 

a 2/3 \ a l/3 J \ a l/3 



W,r',t) = -^-Ai — - Ai -jj t , (28) 



with a = 3eh 2 t 5 / 2 m~ 3 / 2 and Ar", Ar" as above. It follows from Eq. (23) as a short- 
time approximation, to leading order in r. Equation (27), in turn, is recovered if we re- 
place the Airy functions by their asymptotics for large negative argument [64], Ai(— x) — > 

^-1/2(1)1/4^.-3/8 sin (| x 3/2 + zp 

Equation (28) has been obtained combining perturbation theory with semiclassics on 
the level of the van Vleck propagator augmented by a uniform approximation — a route 
that might appear less systematic and controlled than the access through phase-space path 
integration that led to Eq. (23). Its virtue, though, lies in the fact that it is based on 
trajectory pairs, hence readily extends to a higher number of degrees of freedom, while this 
is far from obvious in the case of path integration. This feature together with its simplicity 
and the fact that is is adapted to short propagation time suggests Eq. (28) as a suitable 
choice for on-the-fly molecular- dynamics applications [27, 28]. 

3. Asymptotics e — > 0, K — > 0, t — > 

As pointed out the discussion of Eq. (11), the semiclassical Wigner propagator replaces the 
Liouville propagator by a smooth quantum spot only if the potential is not purely harmonic. 
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FIG. 3: Different versions of the Wigner propagator for a harmonic potential with weak cubic 
anharmonicity as a function of the final phase-space coordinates (q" ,p") at a scaled time cot = 7r/3: 
Wigner propagator in the semiclassical approximation (23) based on the exact quantum result (37) 
(panel a), the van Vleck approximation (27) (b), phase-space path integration (c), and short-time 
version (28) of the same approach (e). Panels d,f are contributions to Eq. (27) corresponding to 
hyperbolic (cosine term, d) and elliptic (sine term, f) trajectory pairs, respectively, cf. Fig. 2b, c. 
Parameter values are fi = 0.01, e = 0.2, m = 1, uo = 1. The underlying classical trajectory (red 
line in all panels) has been launched at (q',p') = (0, —0.1). Color code ranges from red (negative) 
through white (zero) through blue (positive). 

It should converge to the Liouville propagator for vanishing anharmonicity. Indeed, invoking 
the asymptotic expression lim K ^ K~ 1 Ai(x/«:) = S(x) [65] for the Airy functions in Eq. (23), 
we find for a -»■ 0, G w (r",r',t) = 4V35(Ar'!_)S(Ar^). Separating Ap" and Aq" in the 
arguments of the delta functions, this proves to be equivalent to Gw(r" ,r' ,t) = 5(Ar"), the 
classical Liouville propagator. Considering the dependence of a on e, %, and t, this includes 
the classical limit (for subtleties of this limit for Wigner functions and their propagation, see 
[66]) and the limit t — > 0. Moreover, the powers 1, 2, and 2.5, respectively, with which e, h, 
and t appear in a, indicate a hierarchy of the corresponding limits: At finite h, the smooth 
quantum spot already collapses to a delta function in the limit of weak anharmonicity, i.e., 
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for a quantum harmonic oscillator. And even at finite e, taking into account the different 
scaling with time of Aq" and Ap" as they enter Eq. (27), the lateral (£) extension of the 
quantum spot disappears first with t — > while its longitudinal (fj) dimension (in Ap"- 
direction) scales with t as it does with e. 

III. NUMERICAL RESULTS 
A. Models 

1. Morse oscillator 

We choose the Morse oscillator for a detailed study of semiclassical Wigner propagation 
for two reasons: Above all, being prototypical for strongly anharmonic molecular potentials 
and correspondingly complex dynamics [67-70], it is a widely used benchmark for numerical 
methods in this realm [4, 21, 71-73]. Secondly, the Morse oscillator has the remarkable 
advantage, shared only with very few other potentials, that closed analytical expressions are 
available not only for its energy eigenstates [74] but even for their Wigner representations 
[71], a feature that greatly facilitates the comparison with exact quantum- mechanical results 
[43, 71, 75-78] and hence provides a solid basis for an objective test of the performance of 
semiclassical approximations to the propagator. 

The one- dimensional Morse potential [74] 

V Mmsc {q) = D{l-e- aq f (29) 

is determined by the depth D and inverse width a of the potential well. Quantum- 
mechanically, its spectrum is discrete for < E < D and continuous for E > D. The 
number of bound states |a) is given (up to 0(1)) by the parameter A = y/2mD/ah, to 
be interpreted as an inverse Planck's constant in natural units. Analytical expressions in 
terms of Laguerre polynomials are available for the eigenfunctions ip a {<l) = (^l^); see Refs. 
[71, 74]. Applying the Weyl transform (1) to p — \a)(ot\ then leads to the corresponding 
Wigner eigenfunctions W aa (r) [71]. Exact solutions for the continuum states, on the other 
hand, are not known. Numerical evidence indicates, however, that the error caused by their 
omission in the basis set underlying quantum calculations of the time evolution is acceptable 
as long as the energy does not come too close to the threshold E = D. 
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2. Quartic double well 



The quartic double-well potential is in many respects complementary to the Morse os- 
cillator. From a phenomenological point of view, it is the standard model for the study of 
coherent tunneling, and therefore constitutes a particularly hard problem for semiclassical 
propagation. Technically, it is characterized by the absence of a continuum, which facilitates 
quantum calculations. At the same time, no analytical expressions for its eigenfunctions are 
known. 

To begin with, define the quartic double-well potential as V(x) = —mu 2 x 2 /4 + 
m 2 o; 4 x 4 /64£'b, where u is the oscillation frequency near the minima at x± = ±\J&Eb/muo 2 
and E\> their depth. In natural units q = y/mw/hx and r = cut, the Schrodinger equation 
reads 

i|w> = A|*>, ff(M = f-f + lv < 30 > 

Its only parameter, the dimensionless barrier height A = E^/hu, measures approximately 
the number of tunneling doublets, half the number of eigenstates below the barrier top. 

Eigenvalues and eigenfunctions underlying the exact quantum calculations shown in the 
sequel have been obtained in a basis of harmonic-oscillator eigenstates, with unit frequency 
and centered at q — 0. 

B. Propagating delta functions: the propagator as a stand-alone quantity 

To our best knowledge, no detailed account of the propagator of the Wigner function as a 
quantity of its own right has been published to date, except for the recent Ref. [56]. Indeed, 
it is inaccessible to semiclassical approximation schemes based on Gaussian smoothing. We 
therefore find it appropriate to illustrate its basic properties with data obtained for the 
two models featured in this section. Moreover, the interference pattern characterizing the 
Wigner propagator is a sensitive diagnostic for the performance of the different semiclassical 
approximations we are proposing. In addition, an analysis of the propagator allows to 
assess by mere inspection the validity of propagation schemes based on non-classical but 
deterministic trajectories [43-47], see Sec. HID. 

The naked propagator is equivalent to the time evolution of a delta function as initial 
Wigner function, which is not an admissible Hilbert-space element. Notwithstanding, since 



19 



Wigner functions are measurable, e.g., through quantum-state tomography [79], so is the 
propagator, under weak conditions [63], by unfolding the final Wigner function from the 
initial one, which adds legitimacy to considering the propagator as a stand-alone quantity. 
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FIG. 4: Van-Vleck-based semiclassical Wigner propagator (10, 11) (b) and corresponding exact 
quantum result (37) (a) as a function of the final phase-space coordinates (q",p") for the Morse 
oscillator (29) at intermediate energy E = 0.018 at t = 10, approximately the period of the 
underlying classical orbit (red line). Parameters are m = 0.5, D = 1, a = 1.25, fi = 0.005, the 
initial point is (q',p') = (—0.1,0). Color code as in Fig. 3. 

The linear checkerboard structure of the quantum spot appearing in Fig. 3 is owed to 
the cubic anharmonicity of the potential to which the figure refers (its poor resolution in 
the quantum result, panel a, is a numerical artefact). In Fig. 4, we contrast it with data 
for the Morse oscillator, at an intermediate energy E = 0.018.D where the nonlinearity is 
already considerable, and after a propagation time out = 10 of the order of the period of 
the underlying classical orbit. The pattern of Fig. 3 can still be discerned but is severely 
distorted. The nodelines are now strongly curved, and the pattern even folds over, giving 
rise to a phase-space swallow-tail catastrophe. The van-Vleck-based semiclassical propagator 
(10, 11) (panel b) captures this complicated structure surprisingly well. The path-integral- 
based approximation, by contrast, does not reproduce these distortions (not shown). 

Figure 5 shows a similar comparison for the quartic double well, at an energy above the 
barrier top. Here, the interference pattern is even more complex than in Fig. 4. While the 
van-Vleck-based propagator does not agree with the quantum calculation in the intricate 
fine details of the oscillatory fringes, it correctly reproduces much of the more global features 
of the distribution. For an account of its performance at energies below the barrier top, see 
Sec. IIIF. 
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FIG. 5: As Fig. 4 but for the quartic double well (30) with A = 6 at energy E = 2 above the 
barrier at t = 5, approximately the period of the underlying classical orbit (red line). The initial 
point is (q',p') = (0,2). Color code as in Fig. 3. 

C. Propagating Gaussians 

Gaussian initial states deserve their ubiquity not only for their unique physical properties, 
exemplified by coherent states. In the context of the Wigner representation, Gaussians gain 
special relevance as they constitute the only admissible Wigner functions that are positive 
definite and therefore can be interpreted in probabilistic terms [80]. They have achieved a 
fundamental role for semiclassical propagation as they provide a natural smoothing which 
allows to reduce the time evolution of an entire phase-space region to the propagation along 
a single classical trajectory. 

By difference to Gaussian-wavepacket propagation and semiclassical IVRs, however, the 
semiclassical Wigner propagator does not involve any smoothing by construction and thus 
allows to time-evolve arbitrary initial states, Gaussian or not. A more specific question is 
then how semiclassical approximations to it operate on the space of Gaussian phase-space 
distributions. We find that they generally transform them into a broader class of functions 
which, in view of the above theorem, can no longer remain positive definite. 

Define Gaussians in phase space [81] by 



The 2f x 2f covariance matrix A controls size, shape, and orientation of the Gaussian cen- 
tered in r = (po, qo)- The more specific class of minimum-uncertainty Gaussians, equivalent 
to Wigner representations of coherent states [37], is characterized by detA = 1. In what 




(31) 
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follows, in two dimensions r = (p,q), we choose A = diag (1/7,7), so that 



W(r) = — exp 
Tin 



(p-Po) +7 (g-go) 
7/1 



(32) 



equivalent to a position-space wave function [81] ip(q) = (wh) ^ 4 exp[— ^Aq — q \ 2 + ^p 

(q-qo)]- 



1. Autocorrelation functions 



A standard interface between dynamical and spectral data, autocorrelation functions 
have a wide range of applications in atomic and molecular physics and provide a robust 
and easily verifiable assay of the accuracy and efficiency of propagation methods. The 
overlap C^(t) = (ip(t)\ip(0)) for an initial state \ip(0)) upon squaring readily translates 
into an expression in terms of the Wigner representation W^{r) of that state, |C</,(t)| 2 = 
(2ir%y J d 2 ^r W^,(r, 0)W^(r, t), or, involving the propagator explicitly, 



(2iihY J d 2 V j d 2 V^(r",0)G w (r ,/ ,r',t)^(r / ,0). (33) 
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FIG. 6: Autocorrelation function (33) for a Gaussian initial state (31) in a Morse potential (29), 
for the semiclassical approximation (10, 11) (dashed line, red), compared to an exact quantum 
(Eq. (37), full line, blue) and a classical calculation (Eq. (12), dotted, black) of the Wigner prop- 
agator. Parameter values are m = 0.5, D = 1, a = 1.25, h = 0.005, the initial centroid is 
(<7o,Po) = (-0.1,0). 
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FIG. 7: As Fig. 6 but for a quartic double well (30) with A = 6 at an energy E = 1.0026 above 
the barrier top. The initial centroid of the Gaussian is (qo,po) = (4,7). Inset: Nonconservation of 
the norm iV for the same propagation process but without interspersed normalization steps. 

In Fig. 6, we compare the autocorrelation functions (33) for the Morse oscillator with a 
Gaussian initial state (31) as obtained by an exact quantum calculation (full line) to the 
semiclassical approximation (10, 11) (dashed) and the classical propagator (12) (dotted) at 
an intermediate energy. Over the roughly four periods of the classical trajectory monitored, 
the agreement is impressive. Figure 7 is an analogous comparison but for the quartic double- 
well potential. Here, the discrepancy between semiclassical and quantum data is slightly 
more significant. 

In order to obtain the data underlying Figs. 6, 7, it was necessary to periodically renor- 
malize the propagator, compensating for a lack of norm conservation as is notorious for 
van-Vleck-based approximations [82]. In Fig. 7 (inset), we are monitoring the loss of the 
norm of the final Wigner function for a typical run of the propagator (10, 11) without renor- 
malization in the quartic double well over a fraction of the period of the underlying classical 
orbit. 

2. Resolving phase-space structures 

More detailed information on the performance of propagators and reasons for their pos- 
sible failure than is contained in autocorrelation functions can be extracted from the full 
phase-space structures of propagated Gaussians. In particular in strongly nonlinear sys- 
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terns, the degree to which a semiclassical propagator is able to reproduce deviations of the 
time-evolved states from Gaussians is a sensitive measure of its quality. 



r r-r- — i ,j i — — 

;<§)) 

1 1 ' 1 ■ ■' 1 


r 1 — — i ' i 1 — — 

do) 

1 L 1 __l 1 


1 1 1 1 1 1 

- c 


d 


■ e 

: V : 


f 

V : 



-0.2 0.2 -0.2 0.2 

q q 

FIG. 8: Time evolution of a state prepared as a minimum-uncertainty Gaussian (32) (bold blue 
circle in panels a,b is its contour enclosing a Planck cell) in the Morse potential (29) (black contour 
lines in panel a) for the semiclassical approximation (10, 11) (right column), compared to an exact 
quantum calculation (37) (left), at times t = 0.5072 (panels a,b), 1.0144 (c,d), 2.0288 (e,f). The 
full red line is the classical orbit of the initial centroid (qo,po) = (—0.1,0). Parameter values are 
m = 0.5, D = 1, a = 1.25, % = 0.005. Color code as in Fig. 3. 



Such deviations are unmistakable in Fig. 8, a series of snapshots of the time evolution 
of an initially Gaussian state under the same conditions as for the autocorrelation function 
depicted in Fig. 6. The asymmetric droplet shape of the distribution, notably in panels c to 
f, is not compatible with a linearly deformed Gaussian, but is immediately explained by the 
form of the underlying propagator on the centroid trajectory, superimposed on the evolved 
distribution at the same time (Fig. 9). The non-elliptic shape is therefore a consequence of 
the strongly anharmonic potential. In Fig. 8, panels e,f, we even observe fringes where the 
distribution takes on negative values. Unavoidable for non-Gaussian Wigner functions [80], 
this behavior appears even more conspicuous in Fig. 10, analogous to Fig. 8 but for a quartic 
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FIG. 9: Same as Figs. 8c, d but comparing the time-evolved Gaussian (contour lines) to the respec- 
tive propagator (10, 11) on the centroid orbit (red) at the same time t = 1.0144. Color code as in 
Fig. 3. 
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FIG. 10: As Fig. 8 but for the quartic double well (30) with A = 6 at times t = 3 (panels a,b), 9 
(c,d), 15 (e,f). The initial centroid is (qo,Po) = (0,2). Color code as in Fig. 3. 



double well potential. At the same time, the resolution of these structures is relatively poor, 
owing to the reduced density of classical trajectories used in the calculation towards the 
periphery of the initial distribution. 

A more quantitative account of this feature is obtained considering the semiclassical prop- 
agator (23). Folding Airy functions with initial minimum-uncertainty Gaussians leads into 
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a wider function space of polynomials of up to cubic order in the exponent and, in semi- 
classical Wigner propagation, replaces the (frozen, thawed, etc.) Gaussians constituting the 
framework of Gaussian-wavepacket propagation [19-23]. It ranges from strongly distorted, 
oscillatory distributions in the "deep quantum regime" to near Gaussians in the semiclassical 
limit. 

D. Propagating stationary states 

Propagating energy eigenstates of a quantum system may appear a trivial task. For the 
propagation of Wigner functions, however, there is more to gain than just a check of norm 
conservation, as explained in a remarkable work by Lee and Scully [43]. They argue that, 
while for a system with anharmonic potential a phase-space flow along classical trajectories 
cannot reproduce the quantum dynamics, at least in the case of energy eigenstates they 
may be replaced by another type of characteristic: Evidently, Wigner functions associated 
to eigenstates are invariant under a flow that follows their contour lines, suggesting the latter 
as a kind of 'Wigner trajectories" to replace the classical ones for finite h. 

In the present subsection, we confront their proposal with semiclassical Wigner propa- 
gation along pairs of classical trajectories, taking advantage of the fact that also Ref. [43] 
is based on a study of the Morse oscillator. We not only find that our method describes 
the quantum evolution of eigenstates to high accuracy in terms of autocorrelations (far bet- 
ter than [43]) but even present convincing arguments why the very concept of a quantum 
phase-space flow along 'Wigner trajectories" is misleading. 

In Fig. 11 we contrast the contours of the second excited eigenstate of the Morse oscillator 
in Wigner representation [71] with classical trajectories for the same system. At this point, 
it is not even clear which contour to compare with which trajectory. As a reasonable choice 
which becomes compelling in the classical limit, for a given eigenstate |a) we focus on the 
classical trajectory at E — E a . Along this "quantizing" orbit, the corresponding Wigner 
eigenfunctions exhibit a narrow ridge, to be chosen as the relevant contour. We find that 
the classical trajectory roughly follows the contour line but deviates in quantitative detail. 

We would expect semiclassical approximations to improve systematically on mere classical 
propagation. Figure 11 illustrates how this is achieved by propagating along midpoints of 
non-identical classical trajectory pairs: Keeping the initial point r' on a given Wigner contour 
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FIG. 11: Wigner representation of the second excited state of the Morse oscillator [71] (black 
contour lines and color code as in Fig. 3), compared to the classical orbit r^ 2 (t) at the corresponding 
eigenenergy E2 (full red line) and to midpoint paths fj(t) = [rf_(t) + r f (*)]/2 for pairs of classical 
trajectories r^_(f) (dashed red lines), with common initial midpoint f' = r^ 2 (0) but increasing 
initial separation f$ = (p'pq'j) = r'- - r£_ with p'j = and q[ } = 0.002, 0.2, 0.24, 0.26 (see text). 
Parameter values are m = 0.5, D = 1, a = 1.25, h = 0.005. 

fixed, we increase the initial separation f'- = — r'_ of the trajectory pairs launched in 
r'j ± . While for f ' = (classical trajectory) the deviation from the contour line is large, the 
midpoint paths fj(t) = [rj_ (t)+r / - {t)}/2 indeed appear to move towards the Wigner contour 
with increasing f . However, they do not approach it asymptotically but continue shifting 
further past the Wigner contour, indicating that it plays no particular role for quantum time 
evolution in phase space, not even of eigenstates. 

In fact, semiclassical Wigner propagation generally does not reduce to a deterministic 
phase-space flow, not along classical nor along modified (quantum) trajectories. In Fig. 12 
we show the quantum spot formed by the Wigner propagator for the same initial conditions 
as in Fig. 11. As in similar plots throughout this paper, it appears as a smooth distribution 
that cannot be replaced by a delta function, neither on the classical trajectory nor on the 
Wigner contour. It is perhaps surprising but by no means inexplicable that this propaga- 
tor, while reshuffling the Wigner distribution in phase space, nonetheless keeps the Wigner 
eigenfunction invariant. 

This is confirmed by Fig. 13. We are plotting the autocorrelation function for the eigen- 
state Woo( r ) as initial state, which in this case coincides with the norm, propagated with 
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FIG. 12: Van-Vleck based Wigner propagator (10, 11) (panel b) and corresponding exact quantum 
result (37) (a) for the Morse potential (29) as a function of the final phase-space coordinates 
(q",p"), at t = 1.5216. The initial point r' = r^ 2 (0) is the same as in Fig. (11). They are compared 
to the classical orbit (t) (bold line) and the contour of the Wigner eigenfunction W22( r ) (grey), 
cf. Ref. [71], passing through (0). Parameter values are m = 0.5, D = 1, a = 1.25, fi = 0.005. 
Color code as in Fig. 3. 
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FIG. 13: Accuracy of the autocorrelation function for the ground state Wo(r) of the Morse oscillator 
(equivalent to norm conservation) propagated with the semiclassical approximation (10, 11) to the 
Wigner propagator (dashed line, red), a deterministic flow along Wigner contours (full line, green), 
and the classical Liouville propagator (12) (dotted, black). Parameter values are m = 0.5, D = 0.15, 
a = l,h = 1.0. 

the semiclassical Wigner propagator (10, 11), and find deviations of the order of 10~ 5 only, 
as compared to the range of 10~ 4 for propagation along Wigner contours and 10 -3 for mere 
classical Wigner dynamics. Moreover, in Fig. 14 we compare the semiclassically propagated 
state at t = 10 with the initial state. There is no visible difference. We conclude that 




28 



60 
50 
40 
30 

o 20 
g 10 


-10 

-20 

-30 

-0.3 -0.2 -0.1 0.1 0.2 0.3 
P 

FIG. 14: Cross section at p = of the Wigner eigenfunction W22OO) cf. Ref. [71], before (full line, 
blue) and after propagation through a time t = 10 with the semiclassical approximation (10, 11) 
to the Wigner propagator (dashed line, red) and the classical Liouville propagator (12) (dotted, 
black) . 

even in the case of propagating Wigner eigenstates, an acceptable semiclassical approxima- 
tion to the Wigner propagator must not only propagate along non-classical paths but differ 
significantly from a delta function on whatever single trajectory, classical or "quantum" . 

E. Propagating Schrodinger cat states 

The propagation of Schodinger cat states [39] is certainly not a standard task of molecular 
dynamics. Yet it is relevant for the field in various respects: Schrodinger cats are a paradigm 
of quantum coherence and embody the essence of entanglement in a simple setting. They 
allow us to test the performance of propagation methods in this particular respect in an 
objective manner, as the separation of the superposed alternatives and thus the wavelength 
of the corresponding interference pattern can be precisely controlled. In fact, the classic 
double-slit experiment referred to in Ref. [40] to elucidate the challenge quantum coherence 
poses to semiclassical propagation is just a contemporary embodiment of a Schrodinger cat, 
albeit before this term became popular. Finally, a timely, if fancy, application of this notion 
are implementations of quantum computation in medium-size molecules [83]. 

We here consider the propagation of Schrodinger cats prepared as the superposition of 
two coherent states, cf. Eq. (32), with variable separation d initially in position. In the 
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FIG. 15: Schrodinger-cat states time-evolved (34) in the Morse potential (29) at t = 0.3, for 
propagation with the semiclassical approximation (10, 11) to the Wigner propagator (panel b) 
as compared to an exact quantum calculation (37) (a). Parameter values are m = 0.5, D = 1, 
a = 0.25, Ti = 0.005. The initial midpoint and separation, resp., of the Schrodinger cat are 
(^OjPo) = (0-3,0), d = 0.2. Color code as in Fig. 3. 



Wigner representation, this amounts to 



W cat (r) = W(r) + W + (r) + W x (v) 



(34) 



where W±(r) = exp{ — \p\ + 1 2 q±\ h^h} / n , r± = r — [r ± (0, d)], while W x (r) = exp{ — [(p — 
Po) 2 + 1 2 (l — 1o) 2 }/ 1? 1 } x cos[2(p — p )d/h] encodes the quantum coherence in terms of "sub- 
Planckian" oscillations of wavelength %/d'm.p [39]. This initial state is propagated with the 
semiclassical propagator (10, 11) in the Morse potential from the same initial position r as 
in Fig. 8. The result is compared in Fig. 15 to the exact quantum calculation. Apart from 
minor deviations in the shape of the Gaussian envelope, the interference pattern is faithfully 
reproduced. This is not surprising in view of the trajectory-pair construction underlying our 
semiclassical approximation: 

It is instructive to see why propagating along the two classical trajectories of the re- 
spective centroids of the two "classical" Gaussians W±(r) already reproduces essentially 
the sub-Planckian oscillations. The propagator launched from the centroid r of W x (r) 
then comprises two terms, Gw(r",ro,i) = Gwo( r ", r o, t) + G\vx (r", r , t). According to 
Eqs. (10, 11), the first one, propagating along r cl (r ,t), bears no oscillating phase factor 
and therefore practically cancels upon convolution with the strongly oscillatory W x {r'). 
The second one, by contrast, is the contribution of the two centroid orbits r cl (r±,t) form- 
ing a pair of non-identical trajectories. It travels along the non-classical midpoint path 
f x (t) = (r cl (r_, t) + r cl (r + , t))/2 and carries a factor ~ cos[(2/ft)(r + — r_) A r'] = cos[2dp'/h] 
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which couples resonantly to the oscillations in W x {r'). 
F. Tunneling 

Tunneling is to be regarded a quantum coherence effect "of infinite order in H" [84]. One 
therefore does not expect a particularly good performance of semiclassical methods in its 
description, despite various efforts that have been made to improve them in this respect. 
Above all, the complexification of phase space provides a systematic approach to include 
tunneling in a semiclassical framework [5, 85, 86]. Here, by contrast, we restrict ourselves to 
real phase space, in order not to loose the valuable close relationship between Wigner and 
classical dynamics. Even so, we expect that in this framework tunneling can be reproduced 
to a certain degree [87, 88]. To be sure, Wigner dynamics (in real phase space) is exact for 
harmonic potentials. This includes parabolic barriers and hence a specific case of tunnel- 
ing, as pointed out and explained by Balazs and Voros [87]: Since the Wigner propagator 
invariably follows classical trajectories, the explanation rather refers to the initial condition 
in Wigner representation: Owing to quantum uncertainty, it spills over the separatrix even 
if it is concentrated at negative energies, and thus is transported in part with the classical 
flow to the other side of the barrier. 

This is to be considered as a fortunate exception, though. Other, more typical cases 
involving genuine quantum effects, as in particular coherent tunneling between bound states 
as we are considering it here, are not so readily accessible to semiclassical Wigner dynam- 
ics. Quantum tunneling in the Wigner representation, specifically for localized scattering 
potentials, has been studied at depth in [49, 89], however without indicating a promising 
perspective for semiclassical approximations. We are in a slightly more favorite situation as 
the concept of propagation along trajectory pairs provides a viable option how to reproduce 
tunneling by means of a semiclassical Wigner propagator: As illustrated in Fig. 16 (inset), 
it is trajectory pairs with sufficiently separated initial points, probing regions in phase space 
classically inaccessible to one another, which lead to transport along classically forbidden 
paths. 

Figure 16 is the autocorrelation function for a Gaussian initial state prepared at an energy 
slightly below the barrier top of the quartic double well (30). The semiclassical Wigner 
propagator (10, 11) reproduces the revivals but exaggerates their amplitude. This significant 
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FIG. 16: Autocorrelation function (33) for a Gaussian initial state (31) in a quartic double well (30) 
with A = 6 at an energy E = —0.958333 below the barrier top, for the semiclassical approximation 
(10, 11) (dashed line, red), compared to an exact quantum (Eq. (37), full line, blue) and a classical 
calculation (Eq. (12, dotted, black) of the Wigner propagator. The initial centroid of the Gaussian 
is (<7ojPo) = (2,0). Inset: Semiclassical description of coherent tunneling in terms of trajectory 
pairs, in the framework of the van-Vleck based Wigner propagator (10, 11). A wavepacket initially 
prepared near the right minimum of a double-well potential (blue patch) is transported along a 
non-classical midpoint path f(t) = (r_(i) +r^(t))/2 (dashed red line) into the opposite well if the 
two classical orbits r±(t) (full red lines) underlying this path are sufficiently separated initially, 
e.g., on the same side but above the barrier, within the opposite well. Black lines indicate 
other contours of the potential and the separatrix. 

deviation surprises as it overestimates quantum effects. The reasons become clearer upon 
analyzing the full phase-space distribution for the propagated state, Fig. 17, at a time 
corresponding to the right edge of Fig. 16. In terms of global features of the distribution, 
the agreement of the semiclassical (left) with the quantum result (right) is good, remarkably 
even in phase-space regions classically inaccessible from the initial distribution. However, 
while the fine fringes of the quantum distribution reaching out into the opposite well do 
appear in the semiclassical result, their relative weight is not correctly reproduced. This is 
partially explained by the incorrect balance between central peak and oscillatory tail in the 
semiclassical propagator. Moreover, with A = 6, we are in a marginally semiclassical regime 
and cannot expect an optimal performance of the van-Vleck-based Wigner propagator. 
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FIG. 17: Time-evolved state, initially prepared as a minimum-uncertainty Gaussian (32) (bold 
blue circle in panels a,b is its contour enclosing a Planck cell), at time t = 20 in the quartic double 
well (30) with A = 6 at an energy E = —0.958333 below the barrier top, for the semiclassical 
approximation (10, 11) (panel b), compared to an exact quantum calculation (37) (a). The initial 
centroid is (qo,po) = (2,0). Color code as in Fig. 3. 

G. Propagation in the presence of classical chaos 

Besides dynamical coherence effects, complex classical dynamics constitutes a major chal- 
lenge for propagation schemes in molecular physics. The two one-dimensional models dis- 
cussed in the preceding subsections are strongly anharmonic but remain integrable. In order 
to test our method also in the presence of chaos and to check its scalability towards higher 
dimensions, we consider a two-freedom system consisting of Morse oscillators coupled lin- 
early through the positions, a standard model for complex molecular dynamics [90], with a 
potential 

V 2 Morsc(gi, q 2 ) = VMorsc(gi) + ^Morscfe) + Cq^, (35) 

where V^rsefe), i = 1, 2 are Morse potentials (29) and c is the coupling parameter. Starting 
from regular dynamics at c = 0, the system follows the Kolmogorov-Arnol'd-Moser scenario 
[91] with increasing c and becomes fully chaotic for c ^> 1. We here concentrate on the value 
c = 0.3, where phase space is mixed. As initial condition, we choose a two-dimensional 
minimum-uncertainty Gaussian W(r) = W (ri)W (r 2 ) , with W(ri), i = 1,2, as in Eq. (32), 
centered within a major chaotic subregion, see inset of Fig. 18. 

The autocorrelation function depicted in the main part of Fig. 18 includes the first ma- 
jor revival. The van-Vleck-based Wigner propagator reproduces the exact quantum result 
reasonably well and shows a tendency to improve on the classical data, evidencing the role 
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FIG. 18: As Fig. 6 but for two coupled Morse potentials (35) at c = 0.3. Parameter values are 
m = 1, D = 1, a = 1.25, h = 0.125, the initial centroid is (qio,pio, q20,P2o) = (0.5,-0.3,0.4322,0), 
corresponding to a total energy E = 0.5. Inset: Poincare surface of section at P2 = of the 
corresponding classical dynamics. Different initial conditions encoded by colors. The blue circle is 
the contour of the initial Gaussian enclosing a Planck cell. 

of dynamical quantum effects in this system. 

IV. ALGORITHMS 

In view of the objective to demonstrate the viability of semiclassical Wigner propagation 
for numerical applications, we indicate in this section how to construct suitable algorithms 
for this purpose, without entering into details of their implementation. 

A. Assembling the propagator 

1. Exact quantum calculation 

For purposes of comparison, calibration, etc. it is convenient to have an algorithm avail- 
able for a direct exact calculation of the quantum propagator, independent of semiclassical 
approximations. We briefly sketch in the sequel how this is achieved for the propagator of 
the Wigner function. 

Substituting the unitary time-evolution operator U(t) = Y2a e ~ lEat ^ n \ a )( a \ expanded 
in eigenstates \a), H\a) = E a \a), in the definition (4), we obtain the Weyl propagator 
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U w (r,i) = J2 a e lEat / h W a (r) in terms of Wigner eigenf unctions W a (r), Eq. (1) with p = 
|a)(a|. By means of the convolution (5), this yields the Wigner propagator in the form 

(36) 

The convolution in Eq. (36) is an inconvenient feature and can be avoided. Combining 
it with the integrations implicit in the definitions of the Wigner eigenf unctions, it can be 
evaluated, leading to an expression [92] 

G w (r",r',t) = (2^)^e^^-^)'^(r')^(r"), (37) 



which involves generalized Wigner functions [71] 



Pa/3 



q-f ) 



with p a/ 3 = \a){(3\. In the diagonal case a = (3, they correspond to the Wigner representation 
of pure states. For a / ^, they are not generally real as are the diagonal ones (1) but 
Hermitian, W a p{v) = Wg a (r). Both forms, (36) as well as (37), are consistent with the group 
properties of the propagator, in particular with the initial condition Gw(r", r', 0) = 5(r"—r'). 
The quantization required to obtain the basis states |a) can be circumvented by integrating 
the Schrodinger equation directly (e.g., by split-operator methods) to propagate the density 
operator or the wavefunction (for pure states) in a suitable representation, followed by a 
Weyl transform (1). 



2. Semiclassical approximation based on path integrals 

Equations (20, 21) reduce the calculation of the Wigner propagator in the path-integral- 
based approximations to quadratures. More problematic from a practical point of view 
is obtaining the necessary input, since it comprises not only second- but even third-order 
derivatives, see Eq. (21), of the potential — an unavoidable feature, in view of the fundamental 
role the anharmonicity of the potential plays for Wigner propagation. 

This disadvantage is compensated for by the fact that this approximation, at least for 
weak anharmonicity, is accurate enough to allow for a propagation based on a centroid tra- 
jectory only, in analogy to Gaussian wavepacket propagation. It then suffices to evaluate 
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representative values u and a, respectively, of the second and third derivatives of the po- 
tential in the relevant phase-space region, keeping them constant for an entire propagation 
step, Eqs. (23) or (28), till the next update, e.g., in on-the-fly ab-initio molecular dynamics 
[27, 28]. 

3. van-Vleck-based semiclassical approximation 

The Eqs. (10, 11) defining the semiclassical Wigner propagator in van Vleck approxima- 
tion translate into the following straightforward algorithm to compute the propagator as 
such, not operating on any admissible initial Wigner function: 

1. Initial state: Define pairs of initial points r^ ± , j = 1, . . . , N, with common midpoint 
r' = (r'j +r'j _)/2, parameterized, e.g., by spherical coordinates relative to r'. A typical 
value for the number of classical trajectories, used in most of the calculations for one- 
dimensional systems underlying this paper, is N = 10 6 , corresponding to 5 x 10 5 data 
points available for the final coarse-graining, step 3b below. 

2. Time steps: Realize the integration over time as a sequence of L steps tj_i — > ti, 
I = 1,...,L, ti = t' + I At, At = {t" — t')/L. Update the basic ingredients of the 
propagator (10, 11) as follows: 

(a) Trajectories Tj(ti), j = 1, . . . , N, according to the classical force field, 

A rj (t,) = fVHfaitftAt. (39) 

(b) Stability matrices according to the evolution equation M = M^d 2 H(r)/dr 2 (see 
,e.g, [93]), 

AMjft) = M^ ^^ At. (40) 

It suggests itself to implement (a) and (b) as a single step, merging Eqs. (39, 40) 
into a single system of linear equations. 

(c) Actions Sj as (cf. Fig. 1) 

A^(rV) = [r J Jt / )-r J JtO]A[Ar J Jt i )+Ar J+ (tO]/2-[^(r,_(t z ))-i/(r J+ (t z ))]A^ 

3. Final state: 
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(a) Separate elliptic from hyperbolic trajectory pairs according to 

. fell. det [M i+ - M, J > 0, 
traj. pair j is < (41) 

[hyp. det [M j+ - Mj_] < 0. 

(b) Within each of the two sheets, coarse-grain the determinantal prefactor 
|det [M + (t) — M_ (t)] | ~ 1 / 2 and the action S(t) by suitable binning with respect 
to r". In this step, a possibly inhomogeneous distribution of the initial points 
(as, e.g., for polar coordinates) must be accounted for in terms of weight factors. 

(c) Calculate the propagator (10) for each sheet and superpose the two contributions. 



B. Propagating smooth localized initial states: towards Monte Carlo algorithms 

For the more common task of propagating well-localized but quantum-mechanically ad- 
missible initial states (e.g., Gaussians), the method described in the previous subsection is 
not optimal. We can take advantage of the fact that a common midpoint of all trajectory 
pairs is not specified, by evaluating all the N(N — l)/2 pairs formed by a set of N classical 
trajectories to gain a factor O(iV) in efficiency We have to take into account, however, that 
the distribution of centers = (i\, + rfc)/2 of an ensemble of "satellite" phase-space points 
Tj, distributed at random or on an ordered grid with probability density p sat (r), is not this 
density but its self-convolution, p ctr (f) = / d 2 ^r p sat (f — r/2)p sat (r + r/2). For Gaussians 
(31) this reduces to a contraction by a factor y/2. 

Accordingly, we propose the following scheme for the propagation of smooth localized 
initial states: 

1. Initial state: Define a set of initial points r'-, j = 1, . . . , iV sat . A typical value iV sat = 
1000 now generates N ctr = 5 x 10 5 final data points. This can be done in two ways: 

(a) Generate a swarm of random phase-space points r'- covering approximately the 
same phase-space region as the intended initial Wigner function W c t r (r'), and 
associate the weight W ctr (f'j k , t') to each pair with midpoint f'- k , j — 1, . . . , iV sat , 
k = l,...,j-l. 
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FIG. 19: Preparation of initial points for the propagation algorithm for smooth initial distributions, 
Subsection IV B. An ensemble of random "satellite" points r'-, j = l,...,A^ sat (red dots) which 
serve as initial points of iV sa t classical trajectories give rise to N c t T = A r sa t(A r S at — l)/2 midpoint 
paths Yjk{t) starting in the centers f'- k = (r'j + r' k )/2 (black dots) and define the support of the 
propagator at the final time t. The distribution Pctr(r^) (dark grey) of the centers is that of the 
satellites p sa t(ij) (light) contracted by a factor y/2. 

(b) (For Gaussian initial states only) Find the distribution W sat (r' ,t') with covari- 
ance matrix A sat = A ctr /\/2 that entails the intended Wctr(ijfc, t') as its center 
distribution and generate the r'j according to W sat (r',t') (Fig. 19). 

2. Time steps: Propagate classically all iV sat satellites v'- as described in step 2 above. 

3. Final state: Proceed as in step 3 above for every final midpoint f" k . If option lb above 
has been chosen, assign the corresponding weights W ctr (fj- fe , t') to the midpoints f'L in 
the final coarse-graining. 

Being based on ensembles of random phase-space points distributed according to some ini- 
tial density, this propagation scheme readily integrates in Metropolis-type algorithms. This 
suggests itself particularly in the case of high- dimensional spaces where a direct evaluation 
of the phase-space integrals involved would be prohibitive. 
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V. CONCLUSION 



The present work is intended to trigger an interdisciplinary "technology transfer", pro- 
moting the application in molecular physics of a semiclassical phase-space propagation 
scheme that emerged in quantum chaos. We offer a detailed survey of semiclassical Wigner 
propagation, focussing on its implementation and performance as a numerical tool that pro- 
vides a natural initial-value representation. Two complementary approaches to semiclassical 
approximations of the Wigner propagator are considered, one based on the van-Vleck prop- 
agator, including nondiagonal contributions to orbit sums, the other on phase-space path 
integration combined with an expansion of the phase. The former possesses a clear-cut 
classical skeleton in terms of symplectic geometry, is surprisingly accurate in significantly 
nonlinear systems — more so in the reproduction of the notorious fine nodeline structure of 
the Wigner function than in absolute values — and easily generalizes to high- dimensional 
phase spaces. Two known evils of semiclassical propagation, however, return through the 
back door: caustics, now in phase space, and loss of normalization. 

The path-integral approach, by contrast, provides a very precise description at weak 
yet finite anharmonicity. It resolves caustics in terms of Airy functions but tends to fail at 
strong nonlinearity and does not readily extend to higher dimensions. In addition, we devise 
a hybrid based on the same robust structure of trajectory pairs as the van Vleck approach 
but augmented by a uniform approximation to improve the accuracy in the case of close 
stationary points. It appears particularly suitable for heavy-duty numerical applications, 
as in ab-initio molecular-dynamics simulations [27, 28] where integration times between 
updates of the potential are relatively short but the number of freedoms is large and access 
to second- and higher-order derivatives of the potential is exceedingly costly. 

In view of the relatively little that is known to date about semiclassical Wigner prop- 
agation, our results provide sufficient qualitative and quantitative evidence to invalidate 
two popular connotations: The approach readily captures the time evolution of quantum 
coherence effects, including specifically the propagation of Schrodinger-cat states and the 
reproduction of tunneling processes. In these cases, it is crucial that even trajectory pairs 
with large initial separation be taken into account. At the same time, we found that the 
notion of propagating along any kind of "quantum trajectory" is lacking support. While the 
deterministic delta function on the classical trajectory that constitutes the Liouville prop- 
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agator is replaced by a smooth quantum spot, no enhancement of this spot on any other 
deterministic propagation path is observed. 

Various open ends of this work remain to be explored, of which we mention but a few: 

i. Extending semiclassical Wigner propagation to higher dimensions. It is significantly 
facilitated by the transparent geometric underpinning and by viable options to adapt 
it to Monte-Carlo algorithms. 

ii. Semiclassical Wigner propagation in mixed and even fully chaotic systems. Recent 
work in the context of quantum chaos [58] suggests it should work well, but it awaits 
being tested in molecular-physics applications. 

iii. Complexifying phase space. A promising option how to improve on our trajectory- 
based semiclassical approximations, it would resolve caustics and enable a comprehen- 
sive description of tunneling. 

iv. Including incoherent processes like dephasing and dissipation, possibly at finite tem- 
perature. While this requires major modifications of semiclassical approximations 
made for pure states [94], semiclassical Wigner propagation along trajectory pairs read- 
ily extends to systems with Markovian dissipation [31]. Alternatively, the Feynman- 
Vernon influence-functional theory [95, 96] combines well with phase-space path in- 
tegration to achieve high-accuracy semiclassical propagation even in the presence of 
memory effects [97]. 
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